﻿class verlet
{
private:
	caja *box;
	//Paso de tiempo Dt
	double Dt;
	//Dt2
	double Dt2;
	//Dt/2
	double Dth;
	//Dt2/0.5
	double Dt2h;
	//Uso del termostato
	bool termo;
	//Coeficiente de friccion
	double xgamma;
	//sigma = 2*gamma*mi*kB*T
	double sigma;
	//C1 = Dt*(1-1/2*Dt*gamma)
	double C1;
	//C3=0.5*sigma*Dt
	double C3;
	//D1=1/2*Dt*(1-gamma*Dt)
	double D1;
	//D3=gamma*C1*mi
	double D3;
	//D4=sigma*(1-1/2*gamma*Dt)
	double D4;
	//Para el wiener process Wien=sqrt(3*Dt)
	double wien;

public:

	verlet(caja *bx, double xDt, bool xtermo)
	{
		box=bx;
		Dt=xDt;
		Dt2=Dt*Dt;
		Dth=Dt*0.5;
		Dt2h=0.5*Dt*Dt;
		termo=xtermo;
	}

	verlet(caja *bx, double xDt, double sgamma, bool xtermo)
	{
		xgamma=10.0;//sgamma;
		box=bx;
		Dt=xDt;
		Dt2=Dt*Dt;
		Dth=Dt*0.5;
		Dt2h=0.5*Dt*Dt;
		termo=xtermo;
		sigma = sqrt(2.0*xgamma*box->temp0);
		C1 = Dt*(1.0-0.5*Dt*xgamma);
		C3=0.5*sigma*Dt;
		D1=0.5*Dt*(1.0-xgamma*Dt);
		D3=xgamma*C1;
		D4=sigma*(1.0-0.5*xgamma*Dt);
		wien=sqrt(3.0*Dt);

	}

	//Motor LeapFrog con link cell list
	void  LeapFrogLinkCellList(bool calcSv);
	//Motor Verlet clasico con link cell list
	void  verletLinkCellList(bool calcSv);
	//Motor dinamica browniana de Paul y Yoon (PRE 52,2,2076 1995)
	void PaulYoonEngine(bool calcSv);
	//Wiener process
	double wiener(void);
	//Copia literal del motor de Juan Pablo
	void posvel(void);
	//Motor dinamica browniana de Paul y Yoon (PRE 52,2,2076 1995) con linkCell list
	void PaulYoonEngineLinkCell(bool calcSv);
	void LeapFrogSimple(bool calcSv);


};

void  verlet::LeapFrogSimple(bool calcSv)
{
	//si calcSv = true, se calcula el 'shear stress tensor'
	double dv[3];
	double parke=0.;
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz;
	int mc[3], c;

	//Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			box->part[i].pos[j]+=(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt2h;
			//Calculamos la velocidad en t+Dt/2
			if(termo)
			{
				//Empleamos un reescalado de las velocidades como termostato
				box->part[i].vel[j]=(box->part[i].vel[j])*sqrt(box->temp0/box->temp)+Dth*box->part[i].force[j];
			}
			else
			{
				//No utilizamos termostato
				box->part[i].vel[j]+=Dth*box->part[i].force[j];
			}
		}
		box->boundary(i);
	}
	//Calculamos las fuerzas, epotencial y el virial
	interaccionSimple(box, calcSv);

	//Velocidad en t+Dt
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=(box->part[i].force[j])*Dth;
			//Actualizamos la componenete de la energia cinetica
			//parke+=SQR(box->part[i].vel[j]);
			//Calculamos la deriva media de las velocidades de la caja
			dv[j]+=box->part[i].vel[j];
		}

	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	//Eliminamos la deriva de las velocidades de las particulas
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=box->vel[j];
			//Actualizamos la componenete de la energia cinetica
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			//Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
		}
	}

	if(calcSv)
	{
		//Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy;
		box->svxz=svxz;
		box->svyz=svyz;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

void  verlet::LeapFrogLinkCellList(bool calcSv)
{
	//si calcSv = true, se calcula el 'shear stress tensor'
	double dv[3];
	double parke=0.;
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz;
	int mc[3], c;

	//Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			box->part[i].pos[j]+=(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt2h;
			//Calculamos la velocidad en t+Dt/2
			if(termo)
			{
				//Empleamos un reescalado de las velocidades como termostato
				box->part[i].vel[j]=(box->part[i].vel[j])*sqrt(box->temp0/box->temp)+Dth*box->part[i].force[j];
			}
			else
			{
				//No utilizamos termostato
				box->part[i].vel[j]+=Dth*box->part[i].force[j];
			}
		}
		box->boundary(i);
	}
	//Calculamos las fuerzas, epotencial y el virial
	interaccionLinkCellList(box, calcSv);

	//Reset the head of the link-cell list
	for (int i=0;i<box->ncell3;i++)
	{
		box->head[i] = -1;
	}
	//Velocidad en t+Dt
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=(box->part[i].force[j])*Dth;
			//Actualizamos la componenete de la energia cinetica
			//parke+=SQR(box->part[i].vel[j]);
			//Calculamos la deriva media de las velocidades de la caja
			dv[j]+=box->part[i].vel[j];
			//Coordenadas de la nueva celda en la que se encuentra la particula
			mc[j] = int((box->part[i].pos[j]+0.5*box->lbox[0])/box->lgrid);
			box->part[i].celdaxyz[j]=mc[j];
		}
		
		///Indice de la celda donde se encuentra la particula
		c = int(mc[0]*box->ncell2+mc[1]*box->ncell+mc[2]);
		box->part[i].celda=c;
		///Link to the previous occupant (or EMPTY if you're the 1st)
		box->link[i] = box->head[c];
		///The last one goes to the header"""
		box->head[c] = i;
	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	//Eliminamos la deriva de las velocidades de las particulas
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=box->vel[j];
			//Actualizamos la componenete de la energia cinetica
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			//Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
		}
	}

	if(calcSv)
	{
		//Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy;
		box->svxz=svxz;
		box->svyz=svyz;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

void verlet::PaulYoonEngineLinkCell(bool calcSv)
{
	//Velocidad media de las particulas de la caja
	double dv[3];
	double xwien;
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz, svzz;
	double parke=0.;
	int mc[3], c;
	//Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=svzz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			xwien = wiener();
			box->part[i].pos[j]+=C1*(box->part[i].vel[j])+Dt2h*box->part[i].force[j]+C3*xwien;
			xwien = wiener();
			box->part[i].vel[j]+=D1*(box->part[i].force[j])-D3*(box->part[i].vel[j])+D4*xwien;
		}
		//Condiciones de contorno periodicas
		box->boundary(i);
	}
	//Calculamos las fuerzas en la nueva posicion
	interaccionLinkCellList(box, calcSv);

	//Acualizamos la componente de las velocidades que depende de la fuerza en el nuevo paso de tiempo
	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=Dth*(box->part[i].force[j]);
			dv[j]+=box->part[i].vel[j];
		}
	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	//Reset the head of the link-cell list
	for (int i=0;i<box->ncell3;i++)
	{
		box->head[i] = -1;
	}

	//Eliminamos la velocidad de deriva del sistema
	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=box->vel[j];
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			//Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
			svzz+=(box->part[i].vel[2])*(box->part[i].vel[2]);
		}
		///Creamos la link-cell list
		for(int j=0;j<3;j++)
		{
			mc[j] = int((box->part[i].pos[j]+0.5*box->lbox[0])/box->lgrid);
		}
		c = int(mc[0]*box->ncell2+mc[1]*box->ncell+mc[2]);
		box->part[i].celdaxyz[0]=mc[0];
		box->part[i].celdaxyz[1]=mc[1];
		box->part[i].celdaxyz[2]=mc[2];
		box->part[i].celda=c;
		if(c>box->ncell3 || c<0)
		{
			cout<<"error"<<endl;
		}

		///Link to the previous occupant (or EMPTY if you're the 1st)
		box->link[i] = box->head[c];
		///The last one goes to the header"""
		box->head[c] = i;
	}

	if(calcSv)
	{
		//Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy;
		box->svxz=svxz;
		box->svyz=svyz;
		box->svzz=svzz;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

void verlet::PaulYoonEngine(bool calcSv)
{
	///Motor de Paul con calculo de interacciones entre todas las particulas de la caja
	//(sin cutoff). Comprobar para distintos numeros de particulas 

	//Velocidad media de las particulas de la caja
	double dv[3];
	double xwien;
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz, svzz;
	double parke=0.;
	int mc[3], c;
	//Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=svzz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			xwien = wiener();
			box->part[i].pos[j]+=C1*(box->part[i].vel[j])+Dt2h*box->part[i].force[j]+C3*xwien;
			xwien = wiener();
			box->part[i].vel[j]+=D1*(box->part[i].force[j])-D3*(box->part[i].vel[j])+D4*xwien;
		}
		//Condiciones de contorno periodicas
		box->boundary(i);
	}
	//Calculamos las fuerzas en la nueva posicion
	interaccionSimple(box, calcSv);

	//Acualizamos la componente de las velocidades que depende de la fuerza en el nuevo paso de tiempo
	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=Dth*(box->part[i].force[j]);
			dv[j]+=box->part[i].vel[j];
		}
	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	//Eliminamos la velocidad de deriva del sistema
	for (int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=box->vel[j];
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			//Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
			svzz+=(box->part[i].vel[2])*(box->part[i].vel[2]);
		}
	}

	if(calcSv)
	{
		//Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy;
		box->svxz=svxz;
		box->svyz=svyz;
		box->svzz=svzz;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

void  verlet::verletLinkCellList(bool calcSv)
{
	//double posNew[3];
	double dr;
	double dv[3];
	double parke=0.;
	double posR;
	int mc[3], c;
	//Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz;
	//Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	//Calculamos las fuerzas
	interaccionLinkCellList(box,calcSv);

	//Reset the head of the link-cell list
	for (int i=0;i<box->ncell3;i++)
	{
		box->head[i] = -1;
	}

	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			posR=(box->part[i].pos[j])+(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt*Dt/2.0;
			//Condiciones de contorno periodicas
			posR-=box->lbox[i]*anint(posR/box->lbox[i]);
			dr=posR-(box->part[i].posOld[j]);
			dr-=int(dr/(0.5*(box->lbox[j])))*(box->lbox[j]);
			box->part[i].vel[j]=((dr)/(2.0*Dt));
			box->part[i].posOld[j]=box->part[i].pos[j];
			box->part[i].pos[j]=posR;
		}
		box->boundary(i);
	}

	///Uso del termostato?
	if(termo)
	{
		vscale2(box);
	}
	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			dv[j]+=box->part[i].vel[j];
		}
	}

	//Calculamos la velocidad media del sistema
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	for(int i=0;i<box->npart;i++)
	{
		//Eliminamos la velocidad de deriva del sistema
		for (int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=dv[j];
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			//Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
		}
	
		///Creamos la link-cell list
		for(int j=0;j<3;j++)
		{
			mc[j] = int((box->part[i].pos[j]+0.5*box->lbox[0])/box->lgrid);
		}
		c = int(mc[0]*box->ncell2+mc[1]*box->ncell+mc[2]);
		if(c>box->ncell3 || c<0)
		{
			cout<<"error"<<endl;
		}

		///Link to the previous occupant (or EMPTY if you're the 1st)
		box->link[i] = box->head[c];
		///The last one goes to the header"""
		box->head[c] = i;
		
	}
	if(calcSv)
	{
		//Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy;
		box->svxz=svxz;
		box->svyz=svyz;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

double verlet::wiener(void)
{
	double alea = double(rand())/RAND_MAX;
	double dd=0.02236068;//SQRT(0.0005)

	return 1.732050808*dd*int(3.*alea-1.5);
}

void verlet:: posvel(void)
{
	  //Calculamos las fuerzas
	//interaccionLinkCellList(box,false);
	double dgg=1.0-0.5*xgamma*Dt;
	double fpos[3];
	double fx;
	double vxcm,vycm,vzcm;
	double parke=0.0;
	double dvx2, dvy2,dvz2;
	dvx2=dvy2=dvz2=0.0;
	double **fvel;
	double wi;

	fvel = new double*[box->npart];
	for(int i=0;i<box->npart;i++)
	{
		fvel[i] = new double[3];
	}
	//interaccionLinkCellList(box,false);
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			fpos[j]=0.5*Dt*Dt*(box->part[i].force[j]);
			fvel[i][j]=(box->part[i].force[j])*0.5*Dt*(1.0-xgamma*Dt);
			wi = wiener();
			box->part[i].pos[j]+=box->part[i].vel[j]*Dt*dgg+fpos[j]+0.5*sigma*Dt*wi;
			box->part[i].pos[j]-=box->lbox[j]*anint(box->part[i].pos[j]/box->lbox[j]);
		}
	}
	
	interaccionLinkCellList(box,false);
	
	for(int i=0;i<box->npart;i++)
	{
		double tt2=0.;
		for(int j=0;j<3;j++)
		{
			tt2 += (box->part[j].vel[j])*(box->part[j].vel[j]);
			fx=fvel[i][j]+0.5*Dt*(box->part[i].force[j]);
			wi = wiener();
			box->part[i].vel[j]+=fx-xgamma*Dt*dgg*box->part[i].vel[j]+(sigma)*dgg*wi;

		}
	}

	vxcm=0.0;
	vycm=0.0;
	vzcm=0.0;
	for(int j=0;j<box->npart;j++)
	{
		vxcm+=(box->part[j].vel[0]/double(box->npart));
		vycm+=(box->part[j].vel[1]/double(box->npart));
		vzcm+=(box->part[j].vel[2]/double(box->npart));
	}

	
	for(int j=0;j<box->npart;j++)
	{
		box->part[j].vel[0]-=vxcm ;
		box->part[j].vel[1]-=vycm ;
		box->part[j].vel[2]-=vzcm ;
		parke+=SQR(box->part[j].vel[0])+SQR(box->part[j].vel[1])+SQR(box->part[j].vel[2]);
		dvx2+=SQR(box->part[j].vel[0]);
		dvy2+=SQR(box->part[j].vel[1]);
		dvz2+=SQR(box->part[j].vel[2]);
	}
	//cout<<dvx2/box->npart<<" "<<dvy2/box->npart<<" "<<dvz2/box->npart<<endl;
	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);


	for(int i=0;i<box->npart;i++)
	{
		delete fvel[i];
	}
	delete []fvel;
}